方程式:
$$ \frac{Du}{Dt}= -\frac{\partial \phi}{\partial x} + \nu\nabla^2 u \\ \frac{Dv}{Dt}= -\frac{\partial \phi}{\partial y} + \nu\nabla^2 v \\ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0 $$渦度:$\displaystyle \zeta \equiv \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ を用いて、渦度方程式: $$ \frac{D\zeta}{Dt}=\nu\nabla^2\zeta $$ が得られる。
非圧縮より、流線関数: $$ u=-\frac{\partial \psi}{\partial y}, v=\frac{\partial \psi}{\partial x} $$ が定義できる。$\zeta=\nabla^2\psi$なので、 $$ \frac{D\nabla^2\psi}{Dt}=\nu\nabla^4 \psi $$ が得られる($\psi$のみの式)。これをスペクトル法で解く(詳しいことは深く説明できないので石岡さんの本を読むべし)。
using FFTW
using Plots
function main(endstep)
# parameter
x1 = 4pi/5
y1 = 4pi/5
x2 = 6pi/5
y2 = 6pi/5
sigma = pi/10
dt = 1.0/32.0
freq=32
N = 256
K = 170
I = 2*N
x = [2*pi*k/I for k=0:I-1]
#nu = 2.0e-6
nu = 1.0e-4
# Initial condition of vorticity
zeta0(x,y) = exp( (cos(x-x1)+cos(y-y1)-2)/sigma^2 ) + exp( (cos(x-x2)+cos(y-y2)-2)/sigma^2 ) # I.C. function of zeta
zeta = zeta0.(x', x) # Initial condition of vorticity
akl_z = fft(zeta) # Initial condition of vorticity in wavenumber space
# Arrays of wavenumber
k1 = [i for i=0:N]; k2 = [i for i=-N+1:-1]
k_sub = [k1; k2]
k = zeros(I, I)
l = zeros(I, I)
for i = 1:I
l[i, :] = k_sub
k[:, i] = k_sub
end
k2=k.*k
l2=l.*l
kl=k.*l
k2l2=k2+l2
ik2l2=1.0./k2l2
k2_l2=k2-l2
# Initial condition of streamfunction in wavenumber space
akl = akl_z .*ik2l2
akl[1,1]=0.0
# memory allocation
r=copy(akl)
s=copy(akl)
nensei = similar(r)
uhat = similar(r)
vhat = similar(r)
u = similar(r)
v = similar(r)
Ahat = similar(r)
Bhat = similar(r)
Fhat = similar(r)
jikan = similar(r)
k1=similar(r)
k2=similar(r)
k3=similar(r)
k4=similar(r)
# Array for plotting vorticity
Zetas = Array{Array{Float64,2}}(undef, endstep÷freq+1)
num=1
Zetas[num]=Float32.(zeta)
function RKG!(f, h, y, r, s, args...) # Runge-Kutta-Gill
r .= h.*f(y, args...)
y .+= r.*0.5
s .= h.*f(y, args...)
y .+= (-1+sqrt(0.5)).*r .+ (1-sqrt(0.5)).*s
r .= (0.5-sqrt(0.5)).*r .- s
s .= h.*f(y, args...)
y .+= r .+ (1+sqrt(0.5)).*s
r .= -(sqrt(2.0)/3.0*(1.0+sqrt(0.5))).*r .+(2.0/3.0*(-1.0+sqrt(0.5))).*s
s .= h.*f(y, args...)
y .+= r .+ s./6.0
end
function RK4(f, h, x, args...) # 4段4位のRunge-Kutta
k1 .= h .* f(x, args...)
k2 .= h .* f(x + 0.5 * k1, args...)
k3 .= h .* f(x + 0.5 * k2, args...)
k4 .= h .* f(x + k3, args...)
x = x .+ (k1 .+ 2.0 .* k2 .+ 2.0 .* k3 .+ k4) ./ 6.0
return copy(x)
end
function jikanhatten(akl, nu, k, l) # advection and viscousity
nensei .= nu.* k2l2 .* akl
uhat .= -1.0im .* l .* akl
vhat .= 1.0im .* k .* akl
# zxhat = -1.0im .* k .* (k.^2 + l.^2) .* akl
# zyhat = -1.0im .* l .* (k.^2 + l.^2) .* akl
u .= ifft(uhat)
v .= ifft(vhat)
# zx = ifft(zxhat)
# zy = ifft(zyhat)
# F = u.*zx .+ v.*zy
# Fhat = fft(F)
Ahat.=fft(u.*v)
Bhat.=fft(v.*v.-u.*u)
Fhat.=-(k2_l2).*Ahat.-kl.*Bhat
jikan .= Fhat.*ik2l2 .- nensei
jikan[1,1]=0
return jikan
end
# main loop
for itr =1:endstep
akl[K+1:I-K, :] .= 0.0 + 0.0im
akl[:, K+1:I-K] .= 0.0 + 0.0im
# RKG!(jikanhatten, dt, akl, r,s, nu, k, l)
akl = RK4(jikanhatten, dt, akl, nu, k, l)
if itr%freq==0
num+=1
zeta_plot = real.(ifft((k2l2).*akl))
Zetas[num] = Float32.(zeta_plot)
print(itr÷freq, " ")
end
end
return Zetas
end
np=100
endstep=32*np
@time Zetas=main(1)
@time Zetas=main(endstep);
2.169812 seconds (5.38 M allocations: 480.841 MiB, 6.22% gc time, 86.28% compilation time) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 553.345757 seconds (2.13 M allocations: 413.989 GiB, 5.38% gc time, 0.06% compilation time)
553秒(約9分)かかった。多分8分くらいまでは速くできそう。
1.5年前に作ったときは900秒だったので、多少は改善してる感(ぜったいfortranのほうが速いと思う)
N=length(Zetas)
anim = @animate for t=1:N
h=heatmap(Zetas[t],clim=(0,1.0),c=:jet,title="Vortex")
plot(h, size=(700,600))
end
filename="uzu.gif"
gif(anim, filename, fps=10)
[ Info: Saved animation to C:\Users\*****\Desktop\Julia\Vortex\uzu.gif
using HDF5
h5open("zeta_0to100s.h5", "w") do file
for i=1:101
name="zeta"*string(i-1)
write(file, name, Zetas[i])
end
end
コードの詳細は私まで聞いてください。
2023/05/07